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The energetically driven Ehrlich-Schwoebel (ES) barrier had been generally accepted as the pri- 
mary cause of the growth instability in the form of quasi-regular mound-like structures observed on 
the surface of thin film grown via molecular beam epitaxy (MBE) technique. Recently the second 
mechanism of mound formation was proposed in terms of a topologically induced flux of particles 
originating from the line tension of the step edges which form the contour lines around a mound. 
Through large-scale simulations of MBE growth on a variety of crystalline lattice planes using lim- 
ited mobility, solid-on-solid models introduced by Wolf- Villain and Das Sarma-Tamborenea in 2+1 
dimensions, we propose yet another type of topological uphill particle current which is unique to 
some lattice, and has hitherto been overlooked in the literature. Without ES barrier, our simula- 
tions produce spectacular mounds very similar, in some cases, to what have been observed in many 
recent MBE experiments. On a lattice where these currents cease to exist, the surface appears to 
be scale-invariant, statistically rough as predicted by the conventional continuum growth equation. 
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I. INTRODUCTION 



The growth of crystalline thin films via molecular 
beam epitaxy (MBE) has attracted several interests 
experimentally-^ due to ever growing applications in 
many fields, and theoretically 3 -^ due to its rich surface 
structures. Incoming flux of atoms deposit onto a sub- 
strate in the layer-by-layer fashion generally makes the 
film completely free of defects. Excess energy of the 
adatoms allows them to diffuse along the surface away 
from their initial landing positions. These atoms tend 
to minimize their energy by moving towards sites with 
high coordination number such as those along island step 
edges. The process produces an instability in the growth 
morphology leading to the formation of surface struc- 
tures. 

Over the past decades, there have been numerous ex- 
perimental evidences of a pyramid-like mound morphol- 
ogy with a well-defined mound shape and a selective 
sloped— The origin of such a structure has mainly been 
attributed to the presence of Ehrlich-Schwocbel (ES) 
barrier— ~— which, according to Burton-Cabrera-Frank 
theory ; 16 i 17 hinders atoms from moving down a terrace 
resulting in the tendency for them to prefer to flow 
"uphill." Several investigators have confirmed mound 
formation due to ES barrier through their computer 
simulations *£r— More recently another type of destabi- 
lizing ES step-edge current was discovered which occurs 
from the imbalance between the flow of atoms along a 
step edge within the same terrace toward a kink site and 
the flow toward the kink site over a corner . 22 ' 23 The kink 
ES current, unlike the former kind, causes the motion 
within the same terrace, and thus can only occur on a 



surface with spatial dimension higher than one. Until 
the past decade, it was believed that ES barrier was the 
sole cause of mound morphology. One of the authors 24 
was able to obtain a mound-like structure without imple- 
menting the ES barrier. Unlike the ES mechanism which 
are energy-assisted, this mounding instability is proba- 
bilistic and topological in nature. A net current arises 
via the unbalanced between atoms diffusing up and down 
a step edge, hence the name "step-edge diffusion" (SED) 
current. On a simple cubic lattice, this type of current 
only occurs around a corner or a kink site which is sim- 
ilar to the kink ES current. The mounding instability 
through SED current did not emerge spontaneously, and 
was observed only after the use of the so-called "noise 
reduction technique" to suppress the deposition and nu- 
cleation noisei 2 ^— It is unclear that SED mechanism, 
initially studied in a simple cubic system, always occurs, 
and always leads to mounding structure in all crystalline 
lattice structures. Do other topological currents exist in 
other structurally different crystalline lattices? 

Aside from a qualitative description of mound shapes, 
surface morphology is quantitatively measured in terms 
of its roughness as defined by 



W[L,t)= (h 2 (t))-(h(t)Y 



1/2 



(1) 



where (h n (t)) = L~ d '^2, i h'2{t) is the average of atomic 
heights (to the n power) over all lattice sites of dimen- 
sion d and lateral size L. It is believed that W(L,t) 
exhibits a power-law scaling behavior of the form, 

W(L,t)~L a f(£(t)/L), (2) 
where the scaling function f(x) and the lateral dynamical 
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correlation length £(t) are given by 

/(x)-(f' ^J' and Z(t)~t 1/Z (3) 
II, X » 1, 

with a, /3, and z = a/ (3, as suggested by the dynamic 
scaling theory, being the roughness, growth, and dynam- 
ical exponents respectively. In 1+1 dimensions (one spa- 
tial dimension + one time) , These exponents directly as- 
sociate a given discrete growth model to a universality 
class. Despite impressive success in 1+1 dimensions, the 
computational results of several toy models do not con- 
form to the theoretical predictions in 2+1 dimensions. 
The presence of mound morphology, a missing feature 
in one spatial dimension, suggests that a given model 
may belong to a different universality class in a different 
dimension, thus rendering the universality class concept 
futile^I-^ 

In this paper we propose, in addition to ES barrier, 
a competing mechanism for mound formation as a con- 
sequence of probabilistic terrace currents due to the ge- 
ometry of a film's crystalline structure. We begin, in 
Section |nl by describing the helical boundary conditions 
essential for constructing representations of various crys- 
tal structures. Growth simulations are then performed 
on each of these structures according to a set of diffusion 
rules. Calculations of surface roughness and the critical 
exponents are carried out in Section IlIII The roughness 
exponent, in particular, implies the existence of mound 
morphology on a particular crystalline thin film. In Sec- 
tion IIV1 we work out the probabilistic currents of a few 
crystal structures for illustrative purposes. In addition to 
the SED current, a new type of topological current is dis- 
covered which also causes particles to flow in the uphill 
direction. We discuss the mound formation mechanism in 
connection with the underlying continuum growth equa- 
tion in Section|Vl and summarize our work in Section lVTl 

II. MODELS AND CRYSTAL STRUCTURES 

Limited mobility solid-on-solid diffusion models as 
means to model the growth of epitaxial thin film con- 
tinue to attract many interests despite its simplicity and 
decades of intense investigations both computationally 
and analytically. In these models, boxes representing 
atoms are sprinkled down from atop and subsequently 
relax to their final atomic positions without any voids or 
overhangs according to a given set of rules. We examine 
the limit of low substrate temperature so that deposition 
atoms can move at most to one of the nearest neighboring 
sites before coming to rest. Three most prominent dif- 
fusion rules are (i) Wolf- Villain— (WV) where adatom 
moves to maximize its nearest neighbor bonds, (ii) Das 
Sarma-Tamborenea^ii (DT) where adatom moves to in- 
crease its nearest neighbor bonds provided that the cur- 
rent number of bond is less than two, and (iii) Edward- 
Wilkinson 3 - (EW) where adatom moves to the nearest 



neighbor with the minimal local height. Conventionally 
most modelling works are performed on two-dimensional 
rectangular lattices making them only applicable to sim- 
ple cubic materials. At this stage, we have chosen to 
experiment with only WV and DT models and disregard 
EW model because we do not believe that height min- 
imization would lead to any structural formation. Fig- 
ure [1] gives a schematic diagram of WV and DT diffusion 
rules in one spatial dimension. 




DT model 



(b) 

FIG. 1: One dimensional (a) Wolf-Villain model where an 
atom moves to maximize its coordination number and (b) Das 
Sarma-Tamborenea model where an atom moves to increase 
its coordination number given that the current bond count is 
less than two. An atom falls on one of the shaded regions can 
move in one of the directions specified by an arrow. 

To overcome complications in representing any crys- 
talline structures, we adopt the helical boundary condi- 
tions''' 5 which represents any (/-dimensional lattice using 
a one-dimensional chain. As an example, a site on an 
M X N rectangular lattice with coordinates (i, j) is lo- 
cated at the (iN + j) th element of the chain of length 
MN (counting from zero). Its four nearest neighbors 
at (i ± and (i,j ± 1) are mapped to element num- 
ber (z±l) mod MN and (i±N) mod MN respectively. 
Other structures can be constructed in a similar manner. 
The main advantage of the use of helical boundary con- 
ditions over the conventional ones lies in the flexibility 
in representing any n-dimensional structure using one- 
dimensional chain of an arbitrary length (not restricted, 
e.g., to integer x integer in the case of simple cubic crys- 
tal) resulting in the simplicity of the simulation code. 

In this work we simulate the growth of 6 crystal struc- 
tures on 7 planes: simple cubic (SC), body-centered cubic 
(BCC), face-centered cubic (FCC) on (001) and (111) 
planes, simple hexagonal (SH), ideal hexagonal closed 
pack (iHCP) where the ratio c/a — -y/8/3, and hexagonal 
closed pack (HCP) where c/a < -\/8/3. The structures 
whose substrate plane is not listed are understood to be 
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on the (001) plane, or (0001) in the case of the three 
hexagonal lattices. In all simulations we assume that 
atoms in the substrate and the film are of the same type 
so that the film's crystal orientation is the same as that 
of the substrate and no stress of any kind is produced 
along the interface. Unlike other conventional simula- 
tions where time is measured in units of monolayers, here 
it is measured in unit-cell layers (UL) since each struc- 
ture viewed from a different plane may contain a different 
number of atomic layers per unit cell. 

III. CRITICAL EXPONENTS AND THE 
EXISTENCE OF MOUND MORPHOLOGY 

Unlike the generic non-equilibrium surface growth 
where voids and defects are prevalent, the study of ki- 
netic surface roughening of non-equilibrium growth mod- 
els of the solid-on-solid type remains a subject under 
much scrutiny£i^ The former has been generally accepted 
as belonging to the Kardar-Parisi-Zhang universality 
class, 34 originated in the study of Eden cluster growth' 55 
and diffusion-limited aggregation ) 36 ' 37 with the dominant 
growth direction pointing along the surface normal, giv- 
ing rise to a term proportional to ~ (Wh) 2 . Solid-on- 
solid growth models, replicating MBE growth, offer much 
richer behaviors especially when the spatial dimension is 
higher than one. It has been known that a slight change 
in a diffusion rule results in an alteration of the univer- 
sality class. In 1+1 dimensions, WV and EW both be- 
long to EW universality class (with a = 1/2, /? = 1/4, 
z = 2) while DT model belongs to MBE class (with 
a = 1, P = 1/3, z = 3) despite closer diffusion rule 
to WV than is EW models Based on these results, to- 
gether with some symmetry arguments^ it can be shown 
that both diffusion rules obey the following continuum 
growth equation: 

dh - 

— = v 2 V 2 h - v^h + A 13 V • (V/i) 3 + A 22 V 2 (Vfr) 2 + f] 

(4) 

with v 2 = for DT model by symmetry, and v 2 is very 
small compared with and A22 but positive for WV 
model. The term proportional to V ■ (V/1) 3 is often ne- 
glected (by setting A13 to zero) because it generates the 
v 2 S/ 2 h term upon renormalization. Physically speaking 
this term (with A13 > 0) gives a dissipative effect sim- 
ilar to the v 2 V 2 h term, but at a shorter length scale. 
The stochastic nature is captured by the Gaussian noise 
r](x,t) where (r](x,t)) = 0, and 

(r)(x,t)r](x",t')) = 2D5(x~x')S{t-t'). (5) 

Provided the validity of Eq. (|4| with the coefficients 
having the same sign as in 1+1 dimensions, the DT ex- 
ponents in 2+1 dimensions are found to be ag^* =2/3, 

/3g>nt = X / 5i and z cont = 10 / 3j 3 Th(3 W y model j n 

this dimension is predicted to produce logarithmically 
smooth surface with cc^yv = Pwv = ® (1°§))> while the 



dynamical exponent still obeys power law scaling with 
z wv = 2. A large-scale simulation on a SC substrate by 
Das Sarma et ai, however, gives a contradictory result.— 
They reported that the DT model behaved as if it were 
in the EW universality class which suggests that v 2 in 
Eq. Q is no longer zero in this dimension. They also 
observed mound formation in the WV simulations, in- 
stead of logarithmically flat surface, with a wv = 1, 
/3 WV =1/4, and z wv — 4 after making use of the noise 
reduction technique. This implies that WV model is not 
in the EW universality class in 2+1 dimensions. More- 
over, the mounds tend to be of roughly equal size — an 
apparent deviation from being scale invariant. 

Recently Haselwandter et al£&£2- have shown that the 
unstable growth observed in WV model could be ex- 
plained using renormalization group approach. They de- 
rived Eq. ((4]) from a master equation describing the in- 
crement of height at each lattice site according to the 
nearest-neighbor sites. By carefully choosing the regular- 
ization parameter upon taking a continuum limit, they 
were able to obtain the values of the coefficients v 2 , v±, 
A13 and A22- Under repeated RG transformations, these 
values flow differently in d = 1 + 1 and 2 + 1 dimensions. 
In particular the negativity of A13 leads to the change in 
the sign of the diffusion coefficient 1/2 which eventually 
leads to the growth instability in the form of an array of 
islands of lateral size ~ 2tt 2 1 ^4 / ^2 1 - While their anal- 
ysis gives a satisfactory account of the origin of mounds, 
their formulation is still appealed to an atypical regu- 
larization procedure with some dimensional dependency, 
and is not so conveniently extensible to analyze a more 
complicated lattice. In particular it is unclear whether 
the mechanism that gives rise to the growth instability 
is the property of the substrate dimension or of lattice 
geometry. 

To investigate the surface morphologies of the select 
crystal lattices, we perform extensive simulations on 
chains with 100 to 250,000 elements. Due to differences in 
the number of atoms in a unit cell, these numbers trans- 
late to, e.g., the substrate of size 10 x 10 to 500 x 500 in 
the case of SC, and roughly 7 x 7 to 353 x 353 cells in the 
case of BCC(001). In all substrate sizes, the simulations 
are performed until they reach the time step beyond the 
point where the surface roughness W saturates. This 
value ranges from 10 4 UL for 100 elements, up to 10 7 
UL for 250,000 elements. Because of the cross-over be- 
havior of the growth exponent j3 at different timescales, 
its value changes slightly prior to the saturation time. 
The representative value of j3 for a given crystal struc- 
ture is obtained from the asymptotic value of the plot 
between P(L) versus 1/L when L — ► 00. The roughness 
and dynamical exponents a and z are computed from the 
slopes of log(W sa t) and log(i sa t) against the log of lateral 
substrate size L respectively. 

Table U and |TT] show the values of the critical expo- 
nents for DT and WV models respectively. We find that 
the hyper-scaling relation z = a//3 is not generally re- 
spected within the simulation accuracy in both models 
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Das Sarma-Tamborenea model 



structure 



mound 



FCC(lll) 

SH 

iHCP 

FCC(OOl) 

SC 

BCC 

HCP 



0.76 ± 0.04 
0.66 ±0.01 
0.65 ± 0.01 



0.20 ± 0.03 
0.21 ±0.05 
0.22 ± 0.04 



0.63 ±0.02 0.23 ±0.05 

0.62 ± 0.02 0.22 ±0.04 

0.57 ± 0.02 0.24 ±0.05 

0.52 ± 0.02 0.24 ±0.04 



3.3 ±0.2 

3.1 ±0.1 

3.2 ±0.1 
3.2 ±0.1 
3.1 ±0.1 
3.1 ± 0.2 
2.8 ±0.2 



yes(?) 
no 
no 
no 
no 
no 
no 



TABLE I: The critical exponents of the DT model for all 
lattice structures, sorted according to a. 



Wolf- Villain model 



structure 


Q 


P 


z 


mound 


SH 


1.12 ±0.01 


0.25 ± 0.03 


3.9 ±0.3 


yes 


SC 


0.94 ±0.02 


0.23 ± 0.02 


3.6 ±0.1 


yes 


iHCP 


0.86 ±0.02 


0.23 ± 0.02 


3.5 ±0.3 


yes 


FCC(lll) 


0.83 ±0.02 


0.20 ± 0.04 


3.6 ±0.3 


yes 


FCC(001) 


0.57 ±0.03 


0.19 ±0.01 


3.0 ±0.2 


no 


BCC 


0.42 ± 0.03 


0.20 ±0.01 


1.9 ± 0.2 


no 


HCP 


0.37 ±0.03 


0.20 ±0.01 


2.2 ± 0.2 


no 



TABLE II: The critical exponents of the WV model for all 
lattice structures, sorted according to a. 



when a is either too high or too low (such as HCP). 
Albeit some small variations, the values of the growth 
exponent (3 agree with the predicted value (/Sg^' = 1/5) 
from the continuum equation in the case of DT model 
across all lattice structures. The roughness exponent a 
and the dynamical exponent z appear to be slightly less 
than those from the continuum predictions (except for a 
of FCC(lll)). The exponents in the WV case, however, 
do not seem to conform to the theoretical prediction es- 
pecially the dynamical exponent (z^v — 2) which ranges 
from approximately 2 to 4. 

The above discrepancy is removed by noticing the 
last column of Table [TT] which indicates the existence 
of mound-like morphologies on each substrate. In the 
case of BCC and HCP surfaces under WV diffusion rule, 
the surface front appears to be kinetically rough without 
any growth instability. We find a complete agreement 
between the values of the dynamical exponent from the 
simulations and that from the prediction of the contin- 
uum growth equation (z^y = 2). (FCC(001) presents 
an exceptional case. We shall defer its discussion until 
Sec. El) Upon a closer examination of the last column of 
Table HI and HT1 we notice unstable mound-like morpholo- 
gies for those structures with a > 0.66 regardless of the 
diffusion model. (The reason for the question mark in 
the case of FCC(lll) under DT model shall become ev- 
ident at the end of Sec. II V Al and IIVBI ) The separation 
between mound and kinetically rough surfaces at a cer- 
tain value of the roughness exponent has been previously 



observed in the experiment. 40 The fact that the critical 
roughness exponent a c having value greater than 0.66 is 
consistent with the presence of mound morphology from 
the simulations having a large enough value of a of the 
linear or nonlinear continuum growth equations and the 
MBE modellings of a SC lattice without ES barriers.— 

It should be noted that a "rough" surface indicates a 
large value of W(L, t). Since we quantify "roughness" ac- 
cording to Eq. (HJ) , ordered structures such as mounds or 
pyramids, tend to be "rougher" than scale-invariant, ki- 
netically rough surfaces because the mound regions tend 
to be much higher, and the troughs of the hills much 
lower, than the average film height (h(t)}. It is likely that 
a large value of the roughness exponent (a ~ 1) would 
indicate mound morphology on a surface. When a = 1, 
one obtains mounds with slope selection, i.e., mounds 
scale the same way as the lateral substrate dimension. 
For a > 1 (a < 1), mounds tend to grow (shrink) in 
size with a larger substrate. Surfaces that contain visible 
mounds, therefore, have a large value of a, in contrast to 
kinetically rough surfaces of small a which appears flat 
upon taking the thermodynamic limit (L — > oo). 



IV. MOUNDS AND MECHANISM OF MOUND 
FORMATION 

It was well established that both step and kink 
Ehrlich-Schwoebel (ES) barriers could explain the for- 
mation of mound-like structures observed in many MBE 
growth experiments ! 11 ' 41 ^ 2 The former prevents an atom 
on an upper terrace to hop down to a lower terrace, 
while the latter arises from a greater likelihood for an 
atom to move along the edge of a terrace toward a kink 
site than for it to move across a corner to a kink site* 2 - 3 - 
Since then several author a 24 ' 43 have proposed a topolog- 
ically induced probabilistic current known as step-edge 
diffusion (SED) current as an additional cause of mound 
formation. This type of current is physically similar to 
the energetically assisted kink ES current above. Their 
analysis, however, was based on a simple cubic structure. 
It is very unlikely that SED current is the only type of 
probabilistic, topological current in existence. An ar- 
ray of other geometrically more complicated crystalline 
structures could give rise to a new class of geometric, 
probabilistic current. 



A. Lattice-dependent mound morphology 

In our simulations, we see mounds forming since an 
early stage of the growth naturally without any addi- 
tional efforts for all structures with a > 0.66. As time 
progresses, small mounds shift and coalesce into bigger 
ones, similar to what was reported from the WV results 
using the noise reduction technique i 2 ^^ The coarsening 
behavior with mounds of similar sizes during growth im- 
plies that the growth front is not scale invariant. The 
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merging of mounds ends at i sa t where the correlation 
length £(f sa t) is comparable to system's size, and only 
one mound (and one trough) remains. It is interesting 
to note that the growth and coarsening of mounds is not 
stationary; while a large mound subsumes smaller ones 
in order to grow, its tip does not stand still but shifts 
sideways in a series of disappearance and reemergence of 
a peak. 

Figure [2] shows the surface morphologies prior to sat- 
uration times using WV model on chains with 250,000 
elements. This model favors mounds because atoms tend 
to flow toward kink sites which are most likely to have 
the highest coordination numbers. DT model, on the 
other hand, is more inclined to generate a rough sur- 
face because adatoms generally stick to their original 
landings which already have high coordination numbers. 
Even within the same models, mounds do not assume 
the same form. Mounds found on SH and FCC(lll) 
simulations exhibit strong faceted structures. In par- 
ticular, FCC (111) simulations show a striking ensemble 
of triangular pyramids similar to many kinetic Monte 
Carlo (kMC) simulation results at low temperaturesi 44 i 45 
Surfaces of SC and iHCP, on the other hand, only dis- 
play semi-regular hillocks which do not reflect the under- 
lining lattice structure. Three other lattice structures, 
namely FCC(OOl), BCC, and HCP, do not develop visi- 
ble mounds within WV model. 

We do not see the development of mounds in most of 
our DT simulations on the same set of lattices — with a 
notable exception of FCC(lll) as shown in Figure EJa) . 
While the surfaces of other crystalline structures ap- 
pear to be statistically rough, the surface of FCC(lll) 
show round, hemispherical mound morphologies, visu- 
ally very similar to what was observed in DT model with 
ES barrier— or in high-temperature kMC simulations 20 
on a simple cubic lattice. The second highest value of a 
among all lattices simulated using DT model is SH struc- 
ture. Figure EJb) shows the surface of SH prior to sat- 
uration without any apparent mound morphology. The 
surface of other lattices with a smaller value of a exhibits 
a statistically rough interface as traditionally expected 
from a standard DT simulation. On a closer examina- 
tion, however, we find that the development of mounds 
on FCC(lll) surface using DT model is quite different 
from the ones using WV model. Mounds, in this case, do 
not arise from island formations and coarsening of smaller 
mounds into larger ones. Initially the surface appears to 
be statistically rough. As time progresses, the regions 
which are later to become mounds, develop small cracks 
around them. The cracks then deepen, forming narrow 
troughs which gradually enlarge, splitting the original 
surface into many mounds. The perception of mound 
growth is in fact the deepening and broadening of the 
troughs. Finally at late times, small mounds start to 
merge by the progressive disappearance of troughs which 
separate them. We suspect that in this stage, the corre- 
lation length £ (t) dictates the size of each mound (which 
is comparable to the substrate size as the saturation time 



is reached). 



B. Topologically induced uphill currents 

To understand the mechanism of mound formation in 
both models, one should examine the area nearby a ter- 
race edge which separates two flat regions. One com- 
monly accepted explanation as to why an island nu- 
cleation leads to the formation of a large mound-like 
structure is due to the flow of atoms, on average, to- 
wards the mound region resulting in the net "uphill" 
currenti 17 ' 46 " — Without appealing to the use of ES bar- 
rier, we consider a topologically induced uphill current 
in the spirit of SED current. As anticipated, we find 
that all of the lattice structures that develop mounds ap- 
pear to have SED current. To our surprise however, the 
conventional SED current is almost always cancelled by 
local downhill current. We also discover that SH and 
FCC possess yet another type of geometrically induced 
current. Unlike SED current which flows along an edge 
of a terrace towards a kink site, the new current flows 
in the perpendicular direction towards the edge. We be- 
lieve that the reason why this "terrace diffusion" (TD) 
current has never been observed is because in SC, where 
most simulation a 22 ' 43 are based on, an equal and oppo- 
site current flows downhill. The uphill and downhill cur- 
rents thus, on average, cancel each other leaving only 
SED current. It is worth mentioning that TD current is 
analogous to the edge ES current, whereas SED current 
is to kink ES current. The difference is that the edge ES 
current may occur on a one dimensional substrate, while 
TD current is only present on some crystalline lattices in 
two spatial dimensions. 

To illustrate the difference between SED and TD, con- 
sider a step terrace lying along the [1000] direction of a 
SH(001) substrate as shown in Fig. S) Atoms on the up- 
per terrace are denoted by empty circles while those on 
the lower terrace are represented by shaded circles. Ac- 
cording to WV model, a newly deposited atom which falls 
far from the edge of the terrace will not move. The one 
that falls within the proximity of the edge will advance 
along the direction(s) as shown by the arrow(s) in order 
to maximize its bondings. A site with two or more arrows 
indicates that there is an equal probability for an atom 
dropping on it to move in one of the allowed directions. 
Along the flat region away from the kink, atoms tend to 
move uphill as much as they move downhill resulting in 
a net zero flux. Near the kink site, we find that an uphill 
flux tends to occur more often. Note in particular that if 
an atom falls onto position A which situates on the edge, 
it will be attracted toward the kink position B creating a 
small SED current. Since WV diffusion rule only allows 
an atom to move to one of the nearest neighbors, the SED 
current only extends a distance of one atomic position. 
On average, however, a particle does not tend to move 
uphill as a result of this current because there is another 
current flowing downhill in the opposite direction (from 
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FIG. 2: WV simulations of (a) SH, (b) FCC(lll), (c) SC, and (d) iHCP lattices at t = 331132 UL. 



C to B) with the same strength. Nevertheless there is 
a net current in the uphill direction near the corner of 
the terrace edge at position E and F. It is not a SED 
current in the traditional sense since the direction of the 
flow is not along the edge but at an angle towards the 
corner. For the lack of a better word, we shall still refer 
to it as step-edge diffusion current because the current 
still appears in the neighborhood of a kink site and has 
a component parallel to a terrace edge. 

The situation is even more perplexing for a step terrace 
along another high-symmetry direction, the [1200] direc- 



tion, as shown in Fig. [3] Notice that an atom falling onto 
position A next to the kink site does not move toward 
the site. No (traditional) SED current exists along this 
direction. We, however, find a net noncancelling flux of 
currents (thick arrows) flowing perpendicular to the ter- 
race edge in the uphill direction. This flux would serve to 
extend the base of the terrace in a future time step. (We 
shall discuss, in Sec.E] a notable exception of FCC(OOl) 
where a non-zero current does not lead to the formation 
of mounds.) To our knowledge, this type of topological 
current has never been reported in the literature. Near 
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FIG. 3: DT simulations of (a) FCC(lll) and (b) SH at t = 331132 UL. 
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FIG. 4: Local probabilistic currents near a step edge along 
[1000] direction of SH(0001) with a step dent. Shaded circles 
represent atoms on the lower terrace while light circles depict 
those on the upper terrace. An atom dropping on any lattice 
site will move, according to WV rule, to one of sites along 
the corresponding arrow. An atom will not move if it falls 
on a site without an arrow. An atom falling on site A, in 
particular, will be driven towards a kink site B producing 
a small SED current. This current is however cancelled by 
another downhill current from C to A. Global net currents 
are denoted by thick arrows. 



the corner, there also exists a SED current similar to 
those in Fig. |4j Table IIIII gives a summary of the type 
of currents along a given direction during the growth on 
SH, SC, iHCP, FCC(lll), and FCC(001) surfaces. The 
upper terrace resides on the inside of the geometrical fig- 
ures. It is interesting to note that FCC(lll) simulations 
show very strong triangular pyramidal mounds oriented 
in the same direction, and never an inverted triangular 
version. We believe that this is due to the difference be- 




1000J 



FIG. 5: Local probabilistic currents near a step edge along 
[1200] direction of SH(0001). The upper terrace is on the 
left side while the lower terrace on the right. There is a net 
uphill terrace current acting along the [1000] direction. Non- 
cancelling currents are indicated by thick arrows. 



tween the symmetry of the two types of currents; TD 
current is only three-fold symmetric while SED current 
has a six-fold symmetry. The preferred faces are ori- 
ented perpendicular to the directions of the TD currents, 
forming an upright triangular pyramid. The other struc- 
tures whose surface has irregular mounds, namely SC and 
iHCP, are devoid of the TD currents. In addition to the 
edge ES current, TD current should cause an instability 
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forming equilibrium faceting along some vicinal surfaces. 




TABLE III: Edges of SH, SC, iHCP, FCC(lll) and FCC(OOl) 
are shown where the corresponding terrace diffusion (TD) and 
step-edge diffusion (SED) currents are nonzero, calculated 
based on WV diffusion rule. Upper terraces are shown in 
gray. The horizontal direction, indicated by the arrow on the 
upper left hand corner, designates the [fOO] direction for SC 
and FCC(OOf), [fOOO] direction for SH and iHCP, and [ffO] 
direction for FCC(fff). 
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FIG. 6: A step edge along the [100] direction on a BCC(OOI) 
plane separates the lower region (on the bottom of the figure) 
and the upper region (on the top). Three atomic layers are 
present with a lighter shade signifies a higher layer. Occupied 
lattice sites are represented by full circles while dashed circles 
denote unoccupied ones. Locally atoms tend to move away 
from the terrace edge. 

For BCC and HCP, mounds are not observed and nei- 
ther type of current is present. Consider as an example 
the current consideration in the case of BCC(OOl), obey- 
ing WV model. Figure [5] shows three layers of atoms. 
Full circles signify occupied lattice sites. If an atom falls 
on one of these positions, it would have to move along a 
direction designated by one of the arrows towards an un- 
occupied site (dashed circle) . Not only does a net current 
not exist, locally it flows perpendicularly away from the 
terrace edge in both directions. Any atom deposits near 
the edge will likely be pushed away from it which im- 
plies that had an island been formed, its territory would 



not have been extended by this process. In addition, we 
also do not find any current flowing along the edge to- 
wards a kink site. In fact, if an atom falls exactly at the 
kink position, it will diffuse away from the kink. This 
counter-intuitive behavior arises from the fact that, for 
BCC, lattice sites along a terrace edge have the lowest 
coordination numbers. Moving towards these sites in the 
uphill direction from the bottom terrace would reduce 
the number of bonds, in contradiction with the WV dif- 
fusion rule. A closer inspection shows that atoms at the 
bottom of the edge already bond with those at the top. 
An atom which falls on either of these two rows adjacent 
to the edge can only roll away from the edge. The step 
edge in this case serves as a topological barrier preventing 
an atom to cross side. The same situation also happens 
in a HCP lattice whose surface is also mound-free. 

We end this section by giving a brief account of the DT 
simulation results. As discussed at the end of Sec. IIV A[ 
no island formations are observed on the surface of these 
lattices, even in the case of FCC (111) where mounds are 
present. This is consistent with the fact that we do not 
find any non-zero uphill TD or SED currents on any 
structure in any direction. Other than FCC(lll), all 
surfaces appear to be kinetically rough with early-time 
behavior following a power law. 



V. DISCUSSION 

It is clear that the evolution of surface morphology 
depends not only on lattice dimension but also on ma- 
terial's crystal structure. In describing the growth of a 
lattice structure obeying a particular diffusion rule, the 
values (and signs) of the parameters (v2, ^4, A13, A22 and 
D) in the associated continuum growth equation need 
to be adjusted accordingly. The growth morphology are 
primarily categorized into two classes: kinetically rough 
scale-invariant or unstable mounding surface. We find 
that the separation between these two growth regimes 
occurs at the roughness exponent a of around 0.66 re- 
gardless of the prescribed diffusion rule. Further analyt- 
ical study is needed to explain the origin of this magic 
number. For kinetically rough surfaces, the dynamical 
scaling theory seems to give an accurate description of 
the behavior of the growth interface using power laws. 
On the contrary, for destabilized mounding morpholo- 
gies, the growth needs to be described in terms of island 
nucleation and island coarsening. We shall leave the anal- 
ysis of the dynamics of mound coarsening in limited mo- 
bility diffusion models for future worki^i 

Contrary to Ref. [24| we are able to obtain mound mor- 
phology without any noise reduction technique. In our 
simulations we find that mounds are recognizable after its 
lateral size reaches about 100 atomic units. We do not see 
mounds comparable in size to theirs. Our supposition is 
that since we expect the parameters of the correspond- 
ing continuum growth equation to be substrate depen- 
dent, the uphill diffusion term —\v>2\S/ 2 h may overcome 
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the Mullin-type diffusion term — |^4|V 4 /i, which tends to 
suppress small fluctuations, at the length scale given by 
lc ~ VrV^I- This length scale l c , in some crystal 
structure, may be larger than the attempted substrate 
simulation scale, thus, mounds may never be observed. 
In addition our mounds are much more irregular than 
the ones obtained using the noise reduced scheme. Our 
value of a for SC is very close to one which implies that 
mounds have a selective slope in agreement with Ref. [28|. 
We believe that the noise reduction technique, in most 
cases, serves to amplify the mound shape and is not a 
necessary scheme to produce mounds. An exception to 
this observation comes about in the case of FCC(OOl). 
Similar to SH case, we find both SED and TD currents 
acting along a terrace edge — albeit not both in the same 
direction — without any growth instability. We suspect, 
in this case, that either the strength of the uphill current 
is too weak in comparison to the randomness from the 
shot noise, or the substrate size is too limited to see the 
mound formation (L < l c ). More quantitative analysis of 
the current and a thorough investigation of the height- 
height correlation function are necessary to address this 
question. 

A few remarks are in order regarding the topological 
currents. Our observation leads us to believe that the 
mechanism of mound formation within our framework 
is due to both the kink SED current and the straight 
TD current. The latter serves as an extra role in en- 
hancing the faceted structure of mounds on the surface 
where it exists. It is true that other, more complicated, 
edge shapes exist which could cause other geometrical 
currents. In a certain coarse-graining sense, dimples and 
pits can be generated by one kink/corner step similar to 
the ones in Fig. HIE] and |U We s ^ m believe that these cur- 
rents can be largely categorized into curvature-dependent 
versus straight terrace edge type of current. In his re- 
view article^ Krug argued that SED current induces 
Jsed ~ Vk(/i), where k is the local curvature of h(x,t) 
on the plane of the substrate. This results in ~ V 4 /i in 
Eq. . Physically this term emerges as a result of the 
line tension due to the curvature of terrace edge. The 
TD current, on the other hand, appears even when the 
radius of curvature is infinite. Given m = \n ■ Vft| where 
n defines the direction along which the current is active, 
TD current gives rise to the anti-diffusive flow along the 
uphill direction Jtd ~ nm/(l + m 2 ) which is approx- 
imately ~ hm for a small surface sloped Our findings 
suggest that both Jsed and Jtd only act along certain 
preferred directions according to the underlying lattice 
structured (In the case of FCC(lll), for example, the 
TD current may point along one of the following three 
directions: [112], [211] or [121].) 

Our simulations also suggest that, in most structures 
with surface growth instability, local current tends to flow 
towards the bottom of a terrace edge both from the up- 
per and lower terraces. (See, e.g., Fig. @] and [31) On a 
kinetically rough surface such as that of BCC and HCP, 



we observe the local current which flows away from the 
bottom of a terrace edge, similar to what is seen in Fig. El 
Although the effect of this current tends to average out 
on a larger scale, its existence gives rise to the current 
of the form Ji oca i ~ ±n(n • V) 2 /i This translates to a 
new term proportional to (h ■ V) 3 /i in the continuum 
equation. (On a one dimensional substrate, this is sim- 
ply d 3 h/dx 3 .) This term has been previously neglected 
based on the rotation and inversion symmetry about the 
growth direction. Given an anisotropy of each lattice 
structure, we do not believe that the new term should be 
discarded from future investigations. 

Finally we see a different mound formation process on 
FCC(lll) plane under DT diffusion rule. The traditional 
picture of island nucleation, followed by particle accretion 
and mound coarsening may not give an accurate descrip- 
tion of the DT structural formation. Without any net 
uphill currents, we expect a completely different mecha- 
nism at work. Visually, mounds on FCC (111) plane do 
not possess up-down symmetry as the ones obtained us- 
ing WV diffusion rule. This is, however, typical of DT 
growth morphology. We therefore still expect the DT 
terms (— ^V 4 /i and A22 V 2 (Vft.) 2 ) to still be effective in 
the continuum growth equation. In light of the supposed 
anisotropy term that might be present, a complete un- 
derstanding demands a more thorough theoretical inves- 
tigation of the growth equation. 



VI. CONCLUSION 

Through large-scale Monte Carlo simulations, we have 
analyzed MBE growth of thin films on several lattice 
structures based on WV and DT models in 2+1 di- 
mensions. We discovers that at the roughness exponent 
of around 0.66-0.76, the surface morphology of the film 
changes from being kinetically rough with power law scal- 
ing to quasi-regular mound-like structures. Without ES 
barrier, we attribute the morphological difference to the 
appearance of topologically induced, probabilistic parti- 
cle currents. These currents not only arise from the line 
tension along the step edges separating several terraces 
of each mound in the form of SED current, they can also 
emerge perpendicular to flat straight terrace edges in the 
uphill direction in the form of TD current. The latter 
only manifests itself in SH and FCC lattices among sev- 
eral others that we have observed. 
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